work_normalization-etc_rough-draft_wild-type_vary-on-state_antisense

Prepare data for various DGE analyses

Get situated

Set working directory (wd)

Load necessary libraries

Set options

• Use normal numbers instead of default scientific numbers in plots
• Do not limit number of overlaps when including labels in plots

options(scipen=999)
options(ggrepel.max.overlaps = Inf)

Initialize necessary functions

Load in Excel spreadsheet of samples names and variables

The spreadsheet includes Alison’s original sample names; we can use this information to associate the new sample names, which are made up of DESeq2 GLM variables, with the old names, which reflect Alison’s wet-lab, library-prep, etc. work

Load in and process featureCounts tables

Process the sense featureCounts table

Process the antisense featureCounts table

Create tibble: Antisense S. cerevisiae counts, sense K. lactis counts

Record tibble t_fc’s positional information in a GRanges object

pos_info will be used in DESeq2 processing, post-processing, etc.

Perform normalization and run DGE analyses

Perform prep work

Establish table of variables for dds—i.e., a “master” model matrix

• dds stands for “DESeq2 dataset” and is a DESeqDataSet object
• variables for dds are
  • strain
  • state
  • time
  • kit (tcn for “Tecan”, ovn for “Ovation”)
  • transcription (N for “nascent”, SS for “steady state”)
  • auxin
  • timecourse
  • replicate
  • technical

Do prep work for running DESeq2

Make the counts matrix

Make the model matrix

Make the DESeqDataSet, dds

• Use counts_data for the featureCount tallies
• Use col_data for setting up the GLM
• Use pos_info for adding feature metadata, subsequent subsetting, etc.

dds <- DESeq2::DESeqDataSetFromMatrix(
    countData = counts_data,
    colData = col_data,
    design = ~ strain,  # Vary on strain: n3-d vs o-d
    rowRanges = pos_info
)

#  Make a back-up of the DESeqDataSet object
bak.dds <- dds

#  How do things look?
# dds %>% BiocGenerics::counts() %>% head()
# dds@rowRanges
# dds@design
# dds@assays

Prefilter dds

#TODO Let’s keep this in mind and try it if we come to think lowly expressed genes are skewing results.

# threshold <- 1000
# dds_filt <- dds[rowSums(BiocGenerics::counts(dds)) >= threshold, ]
# 
#  Breakdown
#    0 13166
#    1 12822
#    2 12719
#    5 12540
#   10 12358
#   20 12144
#   50 11764
#  100 11403
#  200 10927
#  500 10015
# 1000 8822
# 
# rm(threshold, dds_filt)

I. Run PCA, etc.

Generate non-normalized counts

norm_non <- dds[dds@rowRanges$genome == "S_cerevisiae", ] %>%
    SummarizedExperiment::assay() %>%
    as.data.frame()
norm_non$feature_init <- dds@rowRanges$feature_init[
    dds@rowRanges$genome == "S_cerevisiae"
]

#  Associate non-normalized values with feature metadata
norm_non <- dplyr::full_join(
    norm_non,
    t_fc[(t_fc$genome == "S_cerevisiae"), 1:9],
    by = "feature_init"
) %>%
    dplyr::as_tibble()

Generate normalized counts

Calculate rlog-normalized (unblinded) counts

norm_r <- DESeq2::rlog(
    dds[dds@rowRanges$genome == "S_cerevisiae", ],
    blind = FALSE
) %>%
    SummarizedExperiment::assay() %>%
    as.data.frame()
norm_r$feature_init <- dds@rowRanges$feature_init[
    dds@rowRanges$genome == "S_cerevisiae"
]

#  Associate normalized values with feature metadata
norm_r <- dplyr::full_join(
    norm_r,
    t_fc[t_fc$genome == "S_cerevisiae", 1:9],
    by = "feature_init"
) %>%
    dplyr::as_tibble()

Calculate GeTMM-normalized counts

More details on this relatively new method of normalization, which combines inter- and intra-sample normalization methods and (appears to) perform quite well:

• Baraikdar et al. (Truttmann), Exp Gerontol 2023
• Barrett et al. (Hammarlund), G3 2021
• Bedre, self-published 2023 (most recent update)
• Nelson et al. (Wilkins), Nat Microbiol 2022
• Smid et al. (Sieuwerts), BMC Bioinf 2018
• Walker et al. (Kainer), Comput Struct Biotechnol J 2022
• Zatzman et al. (Shlien), Sci Adv 2022

#  Isolate raw counts for samples of interest
raw <- dds %>%
    SummarizedExperiment::assay() %>%
    as.data.frame()
raw$feature_init <- dds@rowRanges$feature_init

#  Associate non-normalized values with feature metadata
raw <- dplyr::full_join(
    raw,
    t_fc[, c(seq(1,9))],
    by = "feature_init"
) %>%
    dplyr::as_tibble()

#  Calculate counts per kb of gene length (i.e., correct counts for gene
#+ length); gene length is initially in bp and converted to kb
rpk <- ((raw[, 1:5] * 10^3) / raw$length)
rpk[, 6:14] <- raw[, 6:14]

#  Calculate normalization factors using the raw spike-in (K. lactis) counts
norm_KL <- edgeR::calcNormFactors(
    raw[(rpk$genome == "K_lactis"), ][, 1:5]
)

#  Create factor for categories (groups)
model_variables <- stringr::str_split(colnames(rpk[, 1:5]), "_") %>%
    as.data.frame(
        row.names = c(
            "sample", "stage", "day", "kit", "tx", "aux", "tc", "rep", "tech"
        ),
        col.names = paste0("s", c(1:5))
    ) %>%
    t() %>%
    tibble::as_tibble()

group <- factor(
    # Second level is numerator, first level is denominator
    model_variables$sample,
    levels = c("o-d", "n3-d")
)

rm(model_variables)

#  Create edgeR DGEList object composed of S. cerevisiae counts per kb gene
#+ length
dgel <- edgeR::DGEList(
    counts = rpk[rpk$genome == "S_cerevisiae", ][, 1:5],
    group = group
)

#  In the DGEList object, include the normalization factors calculated from
#+ spike-in information
dgel$samples$norm.factors <- norm_KL

#  Check that the normalization factors for each library are appropriately
#+ assigned
dgel$samples

#  Scale the values to counts-per-million
norm_g <- edgeR::cpm(dgel) %>% tibble::as_tibble()
norm_g[, 6:14] <- rpk[rpk$genome == "S_cerevisiae", 6:14]
norm_g

#  Clean up unneeded variables
rm(raw, rpk, norm_KL, group)
rm(dgel)  #TODO Delete dgel? Or use it for trying out DE analyses with edgeR?

Calculate TPM-normalized counts

#  Isolate raw counts for samples of interest
raw <- dds %>%
    SummarizedExperiment::assay() %>%
    as.data.frame()
raw$feature_init <- dds@rowRanges$feature_init

#  Associate non-normalized values with feature metadata
raw <- dplyr::full_join(
    raw,
    t_fc[, 1:9],
    by = "feature_init"
) %>%
    dplyr::as_tibble()

#  Calculate counts per kb of gene length (i.e., correct counts for gene
#+ length or do an "RPK normalization"); then, divide RPK-normalized elements
#+ by the sum of sample RPK divided by one million: this does the actual TPM
#+ normalization
rpk <- tpm <- ((raw[, 1:5] * 10^3) / raw$length)
for (i in 1:ncol(rpk)) {
    tpm[, i] <- (rpk[, i] / sum(rpk[, i] / 1e6))
}

tpm[, 6:14] <- raw[, 6:14]
norm_t <- tpm[tpm$genome == "S_cerevisiae", ]

rm(raw, rpk, tpm)

Run PCA with variously normalized counts

Part 1

#  Make the following code generic --------------------------------------------
#+ ...so that we can try it with different normalization objects (counts
#+ normalized in different ways)
# norm <- norm_non
norm <- norm_r
# norm <- norm_g
# norm <- norm_t


#  Create a PCAtools "pca" S4 object for the normalized counts ----------------
#+ Assign unique row names too
obj_pca <- PCAtools::pca(
    norm[, c(1:5)],
    metadata = dds[dds@rowRanges$genome != "K_lactis", ]@colData
)
rownames(obj_pca$loadings) <- make.names(
    dds[dds@rowRanges$genome != "K_lactis", ]@rowRanges$feature,
    unique = TRUE
)


#  Determine "significant" PCs with Horn's parallel analysis ------------------
#+ See Horn, 1965
horn <- PCAtools::parallelPCA(mat = sapply(norm[, c(1:5)], as.double))

Warning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than availableWarning: more singular values/vectors requested than available

#  Determine "significant" principle components with the elbow method ---------
#+ See Buja and Eyuboglu, 1992
elbow <- PCAtools::findElbowPoint(obj_pca$variance)


#  Evaluate cumulative proportion of explained variance with a scree plot -----
scree <- draw_scree_plot(obj_pca, horn = horn$n, elbow = elbow)
scree

# save_title <- paste0("panel-plot", ".", "scree", ".pdf")
# ggplot2::ggsave(paste0(args$directory_out, "/", save_title), scree)
#TODO Work up some logic for outfile names, location(s) for outfiles, etc.


#  Save component loading vectors in their own data frame ---------------------
loadings <- as.data.frame(obj_pca$loadings)

#  Evaluate the component loading vectors for the number of significant PCs
#+ identified via the elbow method plus two
PCs <- paste0("PC", 1:(as.numeric(elbow) + 2))
top_loadings_all <- lapply(
    PCs, get_top_loadings, x = loadings, z = "all", a = TRUE
)
top_loadings_pos <- lapply(
    PCs, get_top_loadings, x = loadings, z = "pos", a = TRUE
)
top_loadings_neg <- lapply(
    PCs, get_top_loadings, x = loadings, z = "neg", a = TRUE
)

names(top_loadings_all) <-
    names(top_loadings_pos) <-
    names(top_loadings_neg) <-
    PCs
# rm(PCs)
# top_loadings_all$PC1 %>% head(n = 20)
# top_loadings_pos$PC1 %>% head(n = 20)
# top_loadings_neg$PC1 %>% head(n = 20)


#  Analyze positive, negative loadings on axes of biplots ---------------------
#+ Look at the top 15 per axis
images <- list()
mat <- combn(PCs, 2)
for(i in 1:ncol(mat)) {
    # i <- 1
    j <- mat[, i]
    
    PC_x <- x_label <- j[1]
    PC_y <- y_label <- j[2]
    
    images[[paste0("PCAtools.", PC_x, ".v.", PC_y)]] <- plot_biplot(
        pca = obj_pca,
        PC_x = PC_x,
        PC_y = PC_y,
        loadings_show = FALSE,
        loadings_n = 0,
        meta_color = "strain",
        meta_shape = "replicate",
        x_min = -100,
        x_max = 100,
        y_min = -100,
        y_max = 100
    )
    #  x and y ranges for non-normalized counts
    #+ x_min = -200000,
    #+ x_max = 200000,
    #+ y_min = -200000,
    #+ y_max = 200000
    #+
    #+ x and y ranges for rlog-normalized counts
    #+ x_min = -50,
    #+ x_max = 50,
    #+ y_min = -50,
    #+ y_max = 50
    #+
    #+ x and y ranges for GeTMM-normalized counts
    #+ x_min = -50000,
    #+ x_max = 50000,
    #+ y_min = -50000,
    #+ y_max = 50000
    #+
    #+ x and y ranges for TPM-normalized counts
    #+ x_min = -10000,
    #+ x_max = 10000,
    #+ y_min = -10000,
    #+ y_max = 10000
    
    #DECISION For now, go with rlog and TPM; continue to test all four
    
    images[[paste0("KA.", PC_x, ".v.", PC_y)]] <-
        plot_pos_neg_loadings_each_axis(
            df_all = loadings,
            df_pos = top_loadings_pos,
            df_neg = top_loadings_neg,
            PC_x = PC_x,
            PC_y = PC_y,
            row_start = 1,
            row_end = 15,  # 30
            x_min = -0.1,
            x_max = 0.1,
            y_min = -0.1,
            y_max = 0.1,
            x_nudge = 0.02,
            y_nudge = 0.04,
            x_label = x_label,
            y_label = y_label,
            col_line_pos = "black",
            col_line_neg = "red",
            col_seg_pos = "grey",
            col_seg_neg = "grey"
        )
    #  x and y ranges (etc.) for non-normalized counts (messy)
    #+ x_min = -0.75,
    #+ x_max = 0.75,
    #+ y_min = -0.33,
    #+ y_max = 0.33,
    #+ x_nudge = 0.02,
    #+ y_nudge = 0.04,
    #+
    #  x and y ranges (etc.) for rlog-normalized counts (nice and clean)
    #+ x_min = -0.1,
    #+ x_max = 0.1,
    #+ y_min = -0.1,
    #+ y_max = 0.1,
    #+ x_nudge = 0.02,
    #+ y_nudge = 0.04,
    #+
    #+ x and y ranges (etc.) for GeTMM-normalized counts (a bit messy)
    #+ x_min = -0.5,
    #+ x_max = 0.5,
    #+ y_min = -0.5,
    #+ y_max = 0.5,
    #+ x_nudge = 0.04,
    #+ y_nudge = 0.02,
    #+
    #+ x and y ranges (etc.) for TPM-normalized counts (a bit less messy)
    #+ x_min = -0.5,
    #+ x_max = 0.5,
    #+ y_min = -0.5,
    #+ y_max = 0.5,
    #+ x_nudge = 0.04,
    #+ y_nudge = 0.02,
    
    images[[paste0("KA.", PC_x, ".v.", PC_y)]]
}

#  How do things look?
images$PCAtools.PC1.v.PC2

images$KA.PC1.v.PC2$PC_x_pos

images$KA.PC1.v.PC2$PC_x_neg

images$KA.PC1.v.PC2$PC_y_pos

images$KA.PC1.v.PC2$PC_y_neg

# images$PCAtools.PC1.v.PC3
# images$KA.PC1.v.PC3
# images$PCAtools.PC1.v.PC4
# images$KA.PC1.v.PC4
# images$PCAtools.PC2.v.PC3
# images$KA.PC2.v.PC3

Part 2

# for(i in 1:length(names(images))) {
#     # i <- 2
#     vector_names <- names(images) %>% stringr::str_split("\\.")
#     
#     if(vector_names[[i]][1] == "PCAtools") {
#         save_title <- paste0("panel-plot", ".", names(images)[i], ".pdf")
#         ggplot2::ggsave(
#             paste0(args$directory_out, "/", save_title), images[[i]]
#         )
#     } else if(vector_names[[i]][1] == "KA") {
#         save_title <- paste0(
#             "panel-plot", ".", names(images)[i], ".1-x-positive.pdf"
#         )
#         ggplot2::ggsave(
#             paste0(args$directory_out, "/", save_title), images[[i]][[1]]
#         )
#         
#         save_title <- paste0(
#             "panel-plot", ".", names(images)[i], ".2-y-positive.pdf"
#         )
#         ggplot2::ggsave(
#             paste0(args$directory_out, "/", save_title), images[[i]][[2]]
#         )
#         
#         save_title <- paste0(
#             "panel-plot", ".", names(images)[i], ".3-x-negative.pdf"
#         )
#         ggplot2::ggsave(
#             paste0(args$directory_out, "/", save_title), images[[i]][[3]]
#         )
#         
#         save_title <- paste0(
#             "panel-plot", ".", names(images)[i], ".4-y-negative.pdf"
#         )
#         ggplot2::ggsave(
#             paste0(args$directory_out, "/", save_title), images[[i]][[4]]
#         )
#     }
# }
#TODO Work up some logic for outfile names, location(s) for outfiles


#  Plot the top features on an axis of component loading range ----------------
#+ ...to visualize the top variables (features) that drive variance among
#+ principal components of interest
p_loadings <- PCAtools::plotloadings(
    obj_pca,
    components = getComponents(obj_pca, 1),
    # components = getComponents(obj_pca, 1:5),
    rangeRetain = 0.05,
    absolute = FALSE,
    col = c("#785EF075", "#FFFFFF75", "#FE610075"),
    title = "Loadings plot",
    subtitle = "Top 5% of variables (i.e., features)",
    # shapeSizeRange = c(4, 16),
    borderColour = "#000000",
    borderWidth = 0.2,
    gridlines.major = TRUE,
    gridlines.minor = TRUE,
    axisLabSize = 10,
    labSize = 3,  # label_size
    drawConnectors = TRUE,
    widthConnectors = 0.2,
    typeConnectors = 'closed',
    colConnectors = 'black'
) +
    # ggplot2::coord_flip() +
    theme_slick_no_legend
p_loadings

#TODO Work up some logic for saving the plot


#  Evaluate correlations between PCs and model variables ----------------------
#+ Answer, "What is driving biologically significant variance in our data?"
PC_cor <- PCAtools::eigencorplot(
    obj_pca,
    components = PCAtools::getComponents(obj_pca, 1:5),
    metavars = c("strain", "replicate"),
    # metavars = c("strain", "replicate", "sample_name"),
    col = c("#785EF075", "#648FFF75", "#FFFFFF75", "#FFB00075", "#FE610075"),
    scale = FALSE,
    corFUN = "pearson",
    corMultipleTestCorrection = "BH",
    plotRsquared = TRUE,
    colFrame = "#FFFFFF",
    main = bquote(Pearson ~ r^2 ~ correlates),
    # main = "PC Pearson r-squared correlates",
    fontMain = 1,
    titleX = "Principal components",
    fontTitleX = 1,
    fontLabX = 1,
    titleY = "Model variables",
    rotTitleY = 90,
    fontTitleY = 1,
    fontLabY = 1
)

Warning: strain is not numeric - please check the source data as non-numeric variables will be coerced to numericWarning: replicate is not numeric - please check the source data as non-numeric variables will be coerced to numeric

PC_cor

#  Get lists of top loadings for GO analyses ----------------------------------
# for(i in c("PC1", "PC2", "PC3", "PC4")) {
for(i in c("PC1", "PC2")) {
    # i <- "PC1"
    #  Positive
    loadings_pos_PC <- rownames(top_loadings_pos[[i]])[1:500]
    save_title_pos_PC <- paste0(
        "top-500.",
        stringr::str_replace_all(get_name_of_var(loadings_pos_PC), "_", "-"),
        ".", i, ".txt"
    )
    # readr::write_tsv(
    #     dplyr::as_tibble(loadings_pos_PC),
    #     paste0(args$directory_out, "/", save_title_pos_PC),
    #     col_names = FALSE
    # )
    #TODO Work up some logic for location(s) for outfiles
    
    #  Negative
    loadings_neg_PC <- rownames(top_loadings_neg[[i]])[1:500]
    save_title_neg_PC <- paste0(
        "top-500.",
        stringr::str_replace_all(get_name_of_var(loadings_neg_PC), "_", "-"),
        ".", i, ".txt"
    )
    # readr::write_tsv(
    #     dplyr::as_tibble(loadings_neg_PC),
    #     paste0(args$directory_out, "/", save_title_neg_PC),
    #     col_names = FALSE
    # )
    #TODO Work up some logic for location(s) for outfiles
}

II. Run analyses with S. cerevisae features only

Perform size-factor estimation

Here, we subset out the K. lactis features. Thus, we are using only S. cerevisiae features in the size-factor estimation. No control genes are used.

dds_SC <- BiocGenerics::estimateSizeFactors(
    dds[dds@rowRanges$genome != "K_lactis", ]
)
dds_SC@colData

DataFrame with 5 rows and 11 columns
                                          strain    state     time      kit transcription    auxin
                                        <factor> <factor> <factor> <factor>      <factor> <factor>
n3-d_Q_day7_tcn_N_aux-T_tc-F_rep1_tech1     n3-d        Q     day7      tcn             N    aux-T
n3-d_Q_day7_tcn_N_aux-T_tc-F_rep2_tech1     n3-d        Q     day7      tcn             N    aux-T
n3-d_Q_day7_tcn_N_aux-T_tc-F_rep3_tech1     n3-d        Q     day7      tcn             N    aux-T
o-d_Q_day7_tcn_N_aux-T_tc-F_rep1_tech1      o-d         Q     day7      tcn             N    aux-T
o-d_Q_day7_tcn_N_aux-T_tc-F_rep2_tech1      o-d         Q     day7      tcn             N    aux-T
                                        timecourse replicate technical              sample_name
                                          <factor>  <factor>  <factor>                 <factor>
n3-d_Q_day7_tcn_N_aux-T_tc-F_rep1_tech1       tc-F      rep1     tech1 CT8_7716_pIAA_Q_Nascent 
n3-d_Q_day7_tcn_N_aux-T_tc-F_rep2_tech1       tc-F      rep2     tech1 CT10_7718_pIAA_Q_Nascent
n3-d_Q_day7_tcn_N_aux-T_tc-F_rep3_tech1       tc-F      rep3     tech1 CT6_7714_pIAA_Q_Nascent 
o-d_Q_day7_tcn_N_aux-T_tc-F_rep1_tech1        tc-F      rep1     tech1 CT2_6125_pIAA_Q_Nascent 
o-d_Q_day7_tcn_N_aux-T_tc-F_rep2_tech1        tc-F      rep2     tech1 CT4_6126_pIAA_Q_Nascent 
                                        sizeFactor
                                         <numeric>
n3-d_Q_day7_tcn_N_aux-T_tc-F_rep1_tech1   2.350752
n3-d_Q_day7_tcn_N_aux-T_tc-F_rep2_tech1   1.858782
n3-d_Q_day7_tcn_N_aux-T_tc-F_rep3_tech1   2.518477
o-d_Q_day7_tcn_N_aux-T_tc-F_rep1_tech1    0.333973
o-d_Q_day7_tcn_N_aux-T_tc-F_rep2_tech1    0.289205

# n3-d Q N rep1 1.444025
# n3-d Q N rep2 1.119446
# n3-d Q N rep3 1.547506
# o-d Q N rep1 0.772653
# o-d Q N rep2 0.526895

Run DESeq2

Call DESeq2 using default parameters

dds_SC <- DESeq2::DESeq(dds_SC)

#  Check model information
DESeq2::resultsNames(dds_SC)[2]
#  [1] "strain_o.d_vs_n3.d"  #HERE
#+ Thus, the model varies on strain, OsTIR-depletion is the numerator,
#+ Nab3-depletion is the denominator
#+     - Numerator: "top" in MA plots, "right" in volcano plots
#+     - Denominator: "bottom" in MA plots, "left" in volcano plots

Call DESeq2::results()

Make an MA plot that colors features by independent filtering

Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
Please use `linewidth` instead.

Make a volcano plot

III. Run analyses with K. lactis features only

Perform size-factor estimation

Here, we subset out—i.e., remove—the S. cerevisiae features and are thus only analyzing K. lactis features, using all K. lactis features in the size-factor estimation.

dds_KL <- BiocGenerics::estimateSizeFactors(
    dds[dds@rowRanges$genome != "S_cerevisiae", ]
)
dds_KL@colData
# n3-d Q N rep1 0.907373
# n3-d Q N rep2 1.034913
# n3-d Q N rep3 0.944618
# o-d Q N rep1 1.001126
# o-d Q N rep2 1.130992

Run DESeq2

Call DESeq2 using default parameters

dds_KL <- DESeq2::DESeq(dds_KL)

#  Check model information
DESeq2::resultsNames(dds_KL)[2]
#  [1] "strain_o.d_vs_n3.d"  #HERE
#+ Thus, the model varies on strain, OsTIR-depletion is the numerator,
#+ Nab3-depletion is the denominator
#+     - Numerator: "top" in MA plots, "right" in volcano plots
#+     - Denominator: "bottom" in MA plots, "left" in volcano plots

Call DESeq2::results()

Make an MA plot that colors features by independent filtering

Make a volcano plot

IV. Run analyses of S.C. in which all K.L. are controlGenes

Perform size-factor estimation

To run analyses using all K.lactis features as controlGenes, we use a logical vector (a vector composed of elements with values of either TRUE or FALSE) obtained from parsing the rowRanges dataframe within the dds object. In essence, we’re saying, “Return TRUE if the rowRanges variable genome has a value of K_lactis; otherwise, return FALSE.” Then, BiocGenerics::estimateSizeFactors() is using the values associated with those TRUEs to isolate the counts for K. lactis-specific features, in turn using those values to calculate size factors.

#  Recalculate size factors using K. lactis features as `controlGenes`
dds_SC.ctrl_KL <- BiocGenerics::estimateSizeFactors(
    dds,
    controlGenes = (dds@rowRanges$genome == "K_lactis")
)
dds_SC.ctrl_KL@colData
#  Using all of the K. lactis features as `controlGenes`
# n3-d Q N rep1 0.907373
# n3-d Q N rep2 1.034913
# n3-d Q N rep3 0.944618
# o-d Q N rep1 1.001126
# o-d Q N rep2 1.130992

Run DESeq2

Call DESeq2 using default parameters

Since we’ve already calculated the size factors, I think we can exclude K. lactis features from our work from here on out. We have to do some index subsetting to accomplish this (see below).

dds_SC.ctrl_KL <- DESeq2::DESeq(
    dds_SC.ctrl_KL[dds_SC.ctrl_KL@rowRanges$genome != "K_lactis", ]
)

#  Check model information
DESeq2::resultsNames(dds_SC.ctrl_KL)[2]
#  [1] "strain_o.d_vs_n3.d"  #HERE
#+ Thus, the model varies on strain, OsTIR-depletion is the numerator,
#+ Nab3-depletion is the denominator
#+     - Numerator: "top" in MA plots, "right" in volcano plots
#+     - Denominator: "bottom" in MA plots, "left" in volcano plots

Call DESeq2::results()

Make an MA plot that colors features by independent filtering

Make a volcano plot